%% Coefficient estimates
% 2017 2 24
clc; clear all; close all;
load('data_all');
workpath00 = pwd;

addpath('toolbox_xt');
addpath('toolbox_plot');


%% Two panel variables
% xtflag and xtval: combofips
% xtflag2 and xtval2: year
% xtflag3 and xtval3: combofips/year
xtflag3 = xtflag + xtflag2/10000;
xtval3 = unique(xtflag3);

%% Fixed effects estimation - grand time effects
nobs = size(YY,1);
ncity = length(xtval);
nyear = 6;

% temp_y = YY;
% temp_x = [XX, ZZ(:,2), logDD, logDC];
% temp_yd = [year_dummy2, year_dummy3, year_dummy4, year_dummy5, year_dummy6];
% temp_x = [temp_x, temp_yd];
% 
% [d_y, m_y] = xt_demean(temp_y, [ones(size(DD,1),1)], xtflag, xtval);
% d_x = xt_demean(temp_x, [ones(size(DD,1),1)], xtflag, xtval);
% d_b = (d_x'*d_x)\(d_x'*d_y);
% 
% temp_trend1 = d_b(end-4:end);
% 
% [~,m_y2] = xt_mean(m_y, xtflag, xtval);

% figure
% trend1 = repmat(m_y2, 1, nyear) + [zeros(ncity,1),repmat(temp_trend1',ncity,1)];
% hold on
% for i = 1:1:size(trend1,2);
%     plot(xtval2(i),trend1(:,i),'*')
% end
% plot(xtval2, repmat(mean(y),nyear,1)+[0; temp_trend1],'*-k','linewidth',3);
% hold off
% title('Grand time-effects');
% 
% figure
% plot(xtval2, trend1)
% hold on
% plot(xtval2, repmat(mean(y),nyear,1)+[0; temp_trend1],'*-k','linewidth',3);
% hold off

%% Fixed effects - City specific
nobs = size(YY,1);
ncity = length(xtval);
nyear = 6;

% temp_y = YY;
% temp_x = [XX, ZZ(:,2), logDD, logDC];
% temp_yd = [year_dummy2, year_dummy3, year_dummy4, year_dummy5, year_dummy6];
% 
% [d_y, m_y] = xt_demean(temp_y, [ones(size(DD,1),1)], xtflag3, xtval3);
% d_x     = xt_demean(temp_x, [ones(size(DD,1),1)], xtflag3, xtval3);
% d_b = (d_x'*d_x)\(d_x'*d_y);
% dd_y = temp_y-temp_x*d_b + ZZ(:,2)*d_b(15); % add area effects
% 
% [m_y,m_y2] = xt_mean(dd_y, xtflag3, xtval3);
% [~,time_grand] = xt_mean(m_y, xtflag2, xtval2);
% 
% resid_fe = temp_y-temp_x(:,1:14)*d_b(1:14)-temp_x(:,17)*d_b(17);
% 
% m_y3 = xt_mean2d(dd_y, xtflag, xtval, xtflag2, xtval2);
% 
% ddd_y = temp_y-temp_x*d_b;
% [~,time_grand_pure] = xt_mean(ddd_y, xtflag2,xtval2);

%% Random effects with time
nobs = size(YY,1);
xtflag3_str = [num2str(xtflag),num2str(xtflag2)];

%% Load estimation results
% load('results_full_level.mat');

modelind = 6;
islog  = 2;

if islog == 1
%     load('results_full.mat');
elseif islog == 2
%     load('results_full_logp1_A.mat');
load('results_full_logp1_A_s3_REML.mat'); %final model
else
%     load('results_full_level.mat');
end

%% Get coefficients - Fixed part
eval(['temp_est = est_s',num2str(modelind),';']);

% fixed coefficients
% rename nopu                 x1
% rename commercial           x2
% rename industrial           x3
% rename retail               x4
% rename singlefamily         x5
% rename multifamily          x6
% rename office               x7
% rename apartment            x8
% rename holdfordevelopment   x9
% rename holdforinvestment    x10
% rename mixeduse             x11
% rename medical              x12
% rename parking              x13
% rename logsize              x14
[a, b, c] = fixedEffects(temp_est.model);
temp_sel = [1:14, 31, 32 , 23];

tab_fixed = [ ...
  c.Estimate ...
, c.SE ...
, c.tStat ...
, c.pValue ...
];
tab_fixed = tab_fixed(temp_sel,:);

tab_head = {'Name', 'Estimate', 'SE', 't-stat', 'p-value'};
tab_name = {'nopu', 'commercial', 'industrial', 'retail', 'singlefamily', ....
    'multifamily', 'office', 'apartment', 'holdfordevelopment', 'holdforinvestment', ...
    'mixeduse', 'medical', 'parking', 'log size', 'log size^2', 'log size^3', 'log dist. to coast'}';

tab_fixed = [tab_head; [tab_name, num2cell(tab_fixed)]];

%% alpha_jt
tab_head = {'combofips', '2005', '2006', '2007', '2008', '2009', '2010'};
tab_name = num2cell(xtval);

% shrunken
tab_alp_jt = temp_est.alp_jt_mat2;
tab_alp_jt = [tab_head; [tab_name, num2cell(tab_alp_jt)]];

%% Time effects

tab_head = {'combofips', '2005', '2006', '2007', '2008', '2009', '2010'};
tab_time = repmat([0, temp_est.td'], ncity, 1) + temp_est.et;
tab_name = num2cell(xtval);

tab_time = [tab_head; [tab_name, num2cell(tab_time)]]; % first year is always zero

%% Spatial effects

if ~isfield(temp_est, 'alp_j')
    temp_est.alp_j = ones(ncity,1)*nan;
end
if ~isfield(temp_est, 'del_j')
    temp_est.del_j = ones(ncity,1)*nan;
end
if ~isfield(temp_est, 'del2_j')
    temp_est.del2_j = ones(ncity,1)*nan;
end
if ~isfield(temp_est, 'del3_j')
    temp_est.del3_j = ones(ncity,1)*nan;
end
if ~isfield(temp_est, 'del4_j')
    temp_est.del4_j = ones(ncity,1)*nan;
end

tab_sp = [ ...
temp_est.del_j, ...
temp_est.del2_j, ...
temp_est.del3_j, ...
temp_est.del4_j, ...
];

tab_head = {'combofips', 'del_1_j', 'del_2_j', 'del_3_j', 'del_4_j'};
tab_name = num2cell(xtval);
tab_sp = [ tab_head; [tab_name, num2cell(tab_sp)]];

%% Save estimates
workpath = pwd;

% save
% figpath = [workpath, filesep, 'figure_draft', filesep, 'fig01a'];
tabpath = [workpath, filesep, 'figure_draft_rev'];
chk_dir(tabpath);


savefilename = ['tab_parameters_log',num2str(islog),'m',num2str(modelind)];
cd(tabpath);
xlswrite(savefilename, tab_fixed, 'Coeffients', 'A1');
xlswrite(savefilename, tab_sp, 'Spatial function', 'A1');
xlswrite(savefilename, tab_alp_jt, 'alp_jt', 'A1');
xlswrite(savefilename, tab_time, 'Time function', 'A1');
cd(workpath00);





















